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SUMMARY 

The use of the boundary element method for analyzing composite micromechanical behavior is demon- 
strated in this study. Stress-strain, heat conduction, and thermal expansion analyses are conducted using 
the boundary element computer code BEST-CMS, and the results obtained are compared to experimental 
observations, analytical calculations, and finite element analyses. For each of the analysis types, the 
boundary element results agree reasonably well with the results from the other methodologies, with explain- 
able discrepancies. Overall, the boundary element method shows promise in providing an alternative method 
to analyze composite micromechanical behavior. 


INTRODUCTION 

In the analysis of composite materials, the behavior at the micromechanical (constituent) scale is often 
of interest. Several approaches have been used previously in examining micromechanical behavior. These 
include anal ytical methods, in which simplified models are used to yield closed form equations that 
describe effective composite properties, and some micromechanical responses (refs. 1 to 5). Discrete methods 
have also been used previously, in which representative volumes of fiber and matrix are modeled directly 
and analyzed (see ref. 6 as one of many examples). The finite element method (FEM) and boundary ele- 
ment method (BEM) 416 examples of discrete methods. 

The objective of this study is to demonstrate the use of the boundary element method for analyzing 
composite micromechanical behavior as an alternative to other approaches. Specifically, composite micro- 
mechanical analyses have been conducted for several types of problems, and the boundary element results 
have been compared to those obtained from several other methodologies. 


BACKGROUND 

The boundary element method is a discrete analysis method in which discretized integral equations 
are solved over the surface of a domain. This methodology differs from the more familiar finite element 
method in which discretized differential equations are solved over the entire volume of a domain. By only 
requiring the surfaces of a domain to be discretized, the boundary element method can allow a significant 
reduction in the complexity and effort required for modeling. 

To model a composite material, the outer surface of the matrix is discretized. To model the fibers, 
specially formulated “Insert Elements" (ref. 7) are used. For these insert elements, only the centerline of 
the fiber is defined, using nodes and elements, and the fiber radius is defined at each node. Both straight 
and curvilinear fibers are permitted. The fiber surfaces and the variation of the field variables in the plane 
of the fiber cross section are represented analytically within the boundary element formulation. To calcu- 
late the variation of the field variables along the length of the fiber, numerical integration is performed. 

By using the boundary element method to model composite materials, efficiency (from a mesh dis- 
cretization standpoint) can be significantly increased as compared to the finite element method. This 


increase in efficiency is primarily due to the fact that the complexity of discretizing the interior of the 
domain is eliminated. By reducing the number of nodes and elements needed to model a problem, com- 
putational efficiency may also improve. However, due to the fact that commercial finite element codes 
are often optimized for a specific computer system, and the boundary element code that was used for this 
study is a research code, it is not reasonable at this point to compare computational efficiency. 

The BEST-CMS (ref. 7) (Boundary Element Solution Technology-Composite Modeling System) 
computer code was used to conduct the boundary element analyses in this study. The analysis capabili- 
ties available within BEST-CMS include elastostatics, steady state and transient heat conduction, and 
steady state and transient thermoelastics. Several other capabilities of the code are worth mentioning. 

The fiber matrix interface can be idealized as a perfect bond, a linear spring interface, or a nonlinear 
spring-Coulomb-friction interface. The matrix material can be defined to have linear elastic or elastic- 
plastic constitutive behavior. 

There are certain restrictions within the BEST-CMS computer code that should be noted. The analy- 
tical integration in the plane of the fiber cross section assumes that the fiber ends are not free surfaces, 
and thus exposed fiber ends cannot be represented. This restriction results in thermal and mechanical 
loads being transferred through the composite matrix to the end of the fiber. Also, the analytical integra- 
tion of the fiber surface assumes that the fiber has a circular cross section. Particularly for angleplied 
composites, these restrictions result in the intersection of fibers with the matrix outer surface not being 
modeled correctly. 

Finally, rigid plane (multi-point constraint) tying boundary conditions are not specifically available 
within BEST-CMS. To accommodate this limitation, either uniform displacements must be applied to 
nonrestrained faces of a model, or nodal averages of results must be taken over faces of interest in the 
computation of effective composite properties. Another consequence of this limitation is that for angle- 
plied composites appropriate symmetry boundary conditions cannot be applied, and thus various nodal 
constraint boundary conditions must be applied in order to set bounds to the results. 


MATERIAL SYSTEM 

For the analyses in this study, a composite material system composed of SiC (SCS-6) fibers with a 
fiber diameter of 145 /an embedded within a titanium (Ti-15V-3Cr-3Sn-3Al or Ti-15-3 for short) matrix 
was used. The room temperature material properties of the constituents of the material are given in table I 
(ref. 8). A fiber volume fraction of 0.34 was used for the analyses. 

Several important assumptions were applied to the material. Unless otherwise specified, residual 
stresses due to material processing and cooldown were not considered. In addition, a perfect bond was 
assumed between the fiber and matrix. Even though these assumptions do not apply to the actual mate- 
rial, the assumptions allowed the analysis to be simplified. In order to ensure that the assumptions 
remained valid, only linear deformations were considered, and the material was only loaded to a small 
percentage of failure strain. To simplify the thermal analyses, the material properties were assumed not 
to vary with temperature, which also is not the case with the actual material. 


MODEL GENERATION 

To generate the boundary element models for this study, a software interface to the commercial 
preprocessor PATRAN (ref. 9) was developed. This interface, titled COMGEN-BEM (Composite Model 
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Generation-Boundary Element Method) (ref. 10), enables the development of composite boundary element 
models based on user supplied: 


- Fiber volume fraction 

- Fiber diameter 

- Model thickness 

- Mesh density 

- Material properties 

- Fiber orientation angle 

- Loads and boundary conditions 

For this study, one-cell square and four-cell square three-dimensional boundary element models were 
used (fig. 1). The model thickness was set equal to the width. The four-cell square model was shown to 
give more accurate boundary element results, but the one-cell square model was used for cases where 
computational efficiency was more important. Eight noded quadrilateral elements are used to model the 
matrix outer surface, and three noded line elements are used to model the fiber elements. 

It is important to note that the outer surface mesh topology shown in the figures is a mesh of the 
matrix only, not of the fibers. This topology for the matrix mesh was shown to give the most accurate 
results. Also, it should again be emphasized that only the outer surface of the composite matrix needed 
to be modeled. By representing the fibers with line elements, to model angleply composites the fiber line 
elements need only to be rotated. This procedure is much simpler than would be required for finite ele- 
ment modeling, where complex interior meshing would be required. 


STRESS-STRAIN ANALYSES 

The first set of analyses involved studying the stress-strain behavior of the SiC/Ti-15-3 material. 

A four-cell square model was used, and ply orientations of [0] 8 , [90] g , [±30] 2g , and [±45] 2s were considered. 
The original boundary conditions included clamping the back face of the model, applying a uniform tensile 
stress to the front face of the model, and calculating the equivalent strain on the front face. Since the 
front face, as described before, cannot be defined as a rigid plane, the equivalent strain was calculated by 
taking the average of the nodal displacements and dividing that value by the model thickness. The stresses 
were applied up to the approximate point where the material behavior became nonlinear for each ply layup. 
The boundary element results were compared to experimental results obtained by Lerch and Saltsman 
(ref. 8) and, where possible, to results obtained from similar models by the finite element program NASTRAN 
(ref 11). First, the results from the stress-strain analyses for each ply layup are presented, followed by a 
discussion of the results. 

The stress-strain results for the [0] g laminate are shown in figure 2. The BEST-CMS results are 
compared to both finite element and experimental results, and the material was loaded to approximately 
44 percent of its failure strain. For this ply layup, the boundary element results are within 6.5 percent of 
the experimental results at the point of maximum discrepancy, and the finite element results are within 
5 percent of the experimental results. Since the material was found to be linear up to a stress level of 
900 MPa (ref. 8), the slight break seen in the experimental curve is most likely an experimental artifact. 

The stress-strain results for the [90] g laminate, which was loaded to 11 percent of its failure strain, 
are shown in figure 3. The boundary element results are within 1.8 percent of the experimental results, 
while the finite element results are within 3.6 percent of the experimental. The break in the experimental 
curve is most likely due to the start of fiber/matrix debonding (ref. 8), which the analyses do not capture. 
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The stress-strain results for the [±30] 2a and the [±45 ] 2b laminates are shown in figures 4 and 5, 
respectively. The [±30 ] 2b laminate was loaded to 20 percent of its failure strain and the [±45 ] 2b laminate 
was loaded to less than 4 percent of its failure strain. Finite element results were not available for these 
layups, due to the complexity of developing finite element models for angleply laminates. The boundary 
element results for the [±30 ] 2b laminate were within 17.9 percent of the experimental results, and the 
boundary element results for the [±45 ] 2b laminate were within 10.6 percent of the experimental results. 

To attempt to explain the discrepancy between the boundary element results and the experimental 
results, several further studies were conducted on the boundary element models for the [0] 8 and [i30] 2g 
laminates. First, the model was loaded with a cooldown from 700 to 21 °C to account for the residual 
stresses that were applied to the actual material (ref. 8) while maintaining the boundary condition of 
clamping the back face of the model. As can be seen in figures 6 Mid 7, adding residual stresses alone did 
not change the results for the linear range of the stress-strain curve. These results were somewhat expected 
since the experimental results were assumed to have zero stress and zero strain at the beginning of the 
test, thus neglecting any residual stresses. The residual stresses would probably become more significant 
if the fiber-matrix debonding behavior was being examined. 

To examine the effects of the boundary conditions, the nodal constraint boundary condition was 
changed from clamping the back face to applying roller boundary conditions to the back, left and bottom 
faces of the four-cell square model. The cooldown to apply residual stresses was also applied to this problem. 
As can be seen in figures 6 and 7, changing the boundary conditions varied the results slightly for both 
ply orientations, thus serving to provide bounds on the results, but did not change the accuracy of the 
results to any significant extent. 

Finally, the equivalent strain results for the clamped boundary condition with the thermal cooldown 
were examined on an interior face of the model, parallel to the front face, for both ply orientations. For 
this set of examinations, the boundary element results for the [0] 8 laminate improved as compared to 
the results obtained on the end face from within 6.5 percent of the experimental results to within 4.66 percent 
of the experimental results. For the [±30] 2s laminate, the boundary element results improved from within 
17.9 percent of the experimental results to within 11.52 percent of the experimental results. 

The reason for the improvement in the results on the interior faces is most likely due to the BEST- 
CMS assumption that the fiber ends are not free surfaces, thus not allowing exposed fiber ends to be repre- 
sented. The assumption has two primary effects on the results on the end face. First, the loads must be 
transferred through the matrix to the end of the fiber, which can affect the end face results. Second, the 
results on the end face are the response of the matrix material only, while on the interior face the response 
of the fiber is more accurately represented. The fiber end assumption also helps to explain why the boundary 
element results for the [0] 8 laminate were less accurate than those for the [90] 8 laminate. Since the load- 
ing for the [90] 8 laminate is perpendicular to the fibers, the fiber end assumption probably plays a much 
less significant role for this laminate as compared to the [0] 8 laminate. 

The [0] fi and [90] 8 laminate results tended to be more accurate than the angleply laminate results. 
The discrepancy of results is most likely due to the BEST-CMS assumption that the fibers must have a 
circular cross section. The effect of this assumption is that for the angleply laminates, the length of the 
fiber must be cut off sufficiently far away from the outer surface of the matrix to assure a circular cross 
section. The reduction in fiber length causes a reduction in the effective fiber volume fraction, thus 
affecting the results. However, the effect of fiber/matrix debonding still needs to be examined, particu- 
larly for the angleply laminates. Including the debonding behavior will most likely improve the accuracy 
of all of the results, particularly at the higher strain levels, where the material nonlinearities become 

significant. 
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HEAT CONDUCTION ANALYSES 


The second set of analyses involved examining the heat conduction capabilities of BEST-CMS, 
including both steady state and transient heat conduction. For these analyses, the SiC/Ti-15-3 was used 
as the base material. However, the constituent properties were assumed to not vary with temperature. 
Also, for the steady state heat conduction analyses, the thermal conductivity of the fiber was varied so 
that the ratio of fiber to matrix conductivities (K f /K m ) ranged from 2.07 (base material) to 10. 

For all of the heat conduction analyses, a unidirectional ply orientation was used. For the steady state 
problems, a four-cell square model was used. A one-cell square model was used for the transient problem 
in order to reduce the computation time. For the steady state problems, a 100 °C temperature gradient 
was imposed first parallel to the fibers (longitudinal case) and then perpendicular to the fibers (transverse 
case). The equivalent heat flux was then calculated on the low temperature face. The boundary element 
results were then compared to those obtained from NASTRAN finite element results and to results obtained 
from the Composite Cylinders analytical model (ref. 4). While experimental results would have been 
preferable to use as a baseline as opposed to an analytical model, experimental results were not available 
for this problem. For the transient analyses, a 100 °C temperature was applied to the back face of the 
model, and the equivalent temperature on the front face of the model was calculated as a function of time. 

For the steady state heat conduction analyses, the equivalent heat flux is plotted as a function of the 
thermal conductivity ratio for the longitudinal case in figure 8 and for the transverse case in figure 9. 

For the longitudinal case, while the finite element results matched the theoretical results exactly, the 
boundary element results varied from the theoretical by 0.8 percent when the Kf/I^ ratio equaled 2, and 
by 21 percent when the K^/K m ratio equaled 10. 

The cause of the discrepancy of the boundary element results, particularly at the higher conductivity 
values, is most likely related to the variation of the heat flux within each phase (fiber and matrix) that 
each method assumes. The finite element results showed, and the theory assumes, that the flux remains 
constant within each constituent, and is independent of the properties of the other constituents. For 
example, the flux values were constant within the matrix, and did not vary, even as the fiber conduc- 
tivity was changed. The boundary element results, on the other hand, did show that the flux varied within 
each constituent and changed based on the changing of fiber properties. The variation within each phase 
became more pronounced as the fiber conductivity was increased, which accounts for the greater discrep- 
ancy of the results at higher fiber conductivities. Only by studying experimental results, however, can a 
true determination be made as to which method gives the most accurate answers. 

For the transverse steady state heat conduction problem, the boundary element results were within 
7.2 percent of the theoretical results, while the finite element results were within 4 percent of the theoreti- 
cal results. The boundary element results were closer to the theoretical results at higher conductivity 
ratios than was seen in the longitudinal case. This increase in precision was probably due to the fact that 
for the transverse heat conduction, the continuity conditions (temperature and flux must remain constant 
across fiber/matrix boundaries) are the same for both the analytical theory and the boundary element 
method. This set of continuity conditions forces interaction between the constituents to be taken into 
account for all of the analysis methods. 

For the transient heat conduction analyses, the equivalent temperature on the front face of the model 
is plotted as a function of time in figure 10. As can be seen in this figure, the temperature approaches an 
asymptotic value of 100 °C (the original applied temperature) as time is increased. This result is to be 
expected, as with no other boundary conditions applied the model should approach a constant steady 
state temperature. 
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THERMAL EXPANSION ANALYSES 

The final set of analyses involved examining the thermal expansion behavior of the SiC/Ti-15-3 
material system. The material properties of the constituents were assumed not to vary with temperature, 
and the fiber volume fraction was varied from 0.25 to 0.34 to 0.45. A one-cell square model with a unidi- 
rectional ply orientation was used. The back face of the model was clamped, a 100 C temperature load 
was applied to the entire model, and the equivalent strain was calculated on the front face. The bound- 
ary element results were compared to NASTRAN results from equivalent finite element models, and to 
analytical results obtained from a rule of mixtures model (ref. 5). Again, while it would have been prefer- 
able to have experimental results to use as a baseline to compare the computational results against, due 
to the lack of experimental results an analytical model had to be utilized. 

The equivalent strain on the front face is plotted as a function of the fiber volume fraction in fig- 
ure 11. As can be seen from the figure, while the boundary element results and the finite element results 
are within 6.5 percent of each other, there is a significant discrepancy between these results and the theo- 
retical results (around 35 percent for the boundary element models). 

To explain these discrepancies, the case where the fiber volume fraction was 0.34 was further examined 
with additional boundary element analyses. First, the model was lengthened so that the thickness was 
equal to twice the model width. The discrepancy in strains was reduced to 17.45 percent as compared to 
the theory, which indicates that the original model was not long enough to eliminate the effects of apply- 
ing the boundary conditions. By looking at a plane slightly on the interior of the model, the discrepancy 
in strains was reduced to 15.89 percent. The cause of the reduction in the discrepancy is probably due to 
the fact that the fiber ends in BEST-CMS cannot be free surfaces, as discussed before. Finally, the nodal 
constraints were changed for the longer model from clamping the back face to applying rollers to the back, 
bottom and left faces. With these boundary conditions, the strain discrepancy was reduced to 10 percent 
when the end face strains were examined and to 6 percent when the strains on an interior face were exam- 
ined. The roller boundary condition allowed the model to expand freely in all three directions, which 
probably is a more accurate representation of the behavior, which yielded the more accurate results. 


CONCLUSIONS 

The boundary element code BEST-CMS was used to conduct composite micromechanical stress- 
strain, heat conduction, and thermal expansion analyses, with the results being compared to those obtained 
by alternate methodologies. For the stress-strain analyses, the original boundary element results were at 
most within 18 percent of the experimental results. The discrepancies were primarily due to the BEST- 
CMS assumptions that the fibers are not free surfaces and must have circular cross sections. For the heat 
conduction analyses, the boundary element results were found to differ from theoretical and finite element 
results due to BEST-CMS having fewer restrictions on the heat flux variation within a constituent phase. 
For the thermal expansion analyses, the boundary conditions and model length were found to have a sig- 
nificant effect on how close the boundary element results matched theoretical calculations. Overall, the 
specially formulated boundary element methods show promise in providing an alternative approach to 
analyzing composite micromechanical behavior. 

Several areas of future work besides those mentioned already (applying fiber/matrix debonding to 
the stress-strain calculations, etc.) are being undertaken. First, efforts are ongoing to conduct analyses of 
woven composite materials. Due to the complex nature of this composite architecture, the advantages of 
the boundary element method model generation become even more significant in analyzing these materials. 
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Also, elastic-plastic matrix behavior will be incorporated into the stress-strain models, in order to more 
accurately capture the material behavior at high strain levels. 
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TABLE I.— MATERIAL PROPERTIES FOR 


COMPOSITE CONSTITUENT MATERIALS 



SiC 

Ti-15-3 

Density, kg/m 3 
Young’s modulus, GPa 
Poisson’s ratio 
Conductivity, W / mK 
Heat capacity, J/kg K 
CTE, 1/*C 

2989.34 

393 

0.19 

16.74 

1256.10 

2.2x10“* 

4760.81 

88 

0.32 

8.10 

502.44 

8.1x10“* 
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Figure 8.— Effect of fiber conductivity on longitudinal steady state heat flux, [0] 8 laminate, FVR = 0.34, temperature gradient = 
100 deg. C, base material = SiC/Ti-15-3. 
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jre 9. — Effect of fiber conductivity on transverse 
)0 deg. C, base material = SiC/Ti-1 5-3. 


steady state heat flux, [0]$ laminate, FVR = 0.34, temperature gradient = 



Figure 1 0. Variation of equivalent temperature with time, [0]g laminate, FVR = G.34, applied temp. = 100 deg. C, base material - 

SiC/TM5-3. 
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Fiber Volume Ratio 

Figure 1 1 . Effect of FVR on equivalent longitudinal strain, [0]g laminate, temp. Increase = 100 deg. C, base material 

SiC/Tl - 1 5 - 3. 
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